Abstract

Evolve a two-components system and compare the results with the equation (7) in http://link.aps.org/abstract/PRA/v59/pR31. This notebook is meant to be a tool to study pseudo-spin oscillation.

Set up

In [34]:
%matplotlib inline
from __future__ import division # to carry out normal division with just one '/'.
import numpy as np
import trottersuzuki as ts
from matplotlib import pyplot as plt

import plotly.plotly as py
import plotly.graph_objs as go
from plotly.offline import download_plotlyjs, init_notebook_mode, iplot
init_notebook_mode()


dim = 100      # square lattice
length = 25.
delta_x = delta_y = length / dim
periods = [0, 0]      
kernel_type = "cpu"      #the evolution is performed by the CPU kernel

# use harmonic oscillator units
particle_mass_A = particle_mass_B = 1.
w_x = w_y = 1.   # frequencies of the harmonic confinement
rot_coord_x = rot_coord_y = dim / 2   # the axis of rotation will be located at the lattice center

def harmonic_pot(x,y):
    return 0.5 * particle_mass_A * (w_x**2 * x**2 + w_y**2 * y**2)

TF_radius = 5.938   # wished TF radius, in harmonic oscillator units
coupling_const = np.pi * (TF_radius**4) / 4.   # this will be multiplying a wavefunction normalized to unity

# vectors containing the physical coordinates (in h.o. units, and centered in the lattice)
x_vec = y_vec = (np.arange(dim) - dim*0.5) * delta_x

# define the external potential
external_pot_A = np.zeros((dim, dim))
for i in range(0, dim):
    y = y_vec[i]
    for j in range(0, dim):
        x = x_vec[j]
        external_pot_A[i,j] = harmonic_pot(x,y)
        
external_pot_B = external_pot_A



# plot energy and norm
def plot_evolution():
    plt.plot(Time, Tot_energy, label = 'Total Energy')
    plt.legend()
    plt.xlabel('Time')
    plt.show()
    plt.plot(Time, Norm_A, label = 'Squared norm A')
    plt.plot(Time, Norm_B, label = 'Squared norm A')
    plt.plot(Time, Norm_A + Norm_B, label = 'Total squared norm')
    plt.legend()
    plt.xlabel('Time')
    plt.axis([0, max_it * iterations * delta_t, 0, np.amax(Norm_A + Norm_B) + 1])
    plt.show()

# exports to file the energies at a given time
def get_info(time):
    tot_energy = ts.calculate_total_energy_2GPE(p_real_A, p_imag_A, p_real_B, p_imag_B, particle_mass_A, particle_mass_B,
                                                [g_A, g_B, g_AB, Omega_Rabi_real, Omega_Rabi_imag], external_pot_A, external_pot_B,
                                                omega, rot_coord_x, rot_coord_y, delta_x, delta_y)
    norm_A = ts.calculate_norm2(p_real_A, p_imag_A, delta_x, delta_y)
    norm_B = ts.calculate_norm2(p_real_B, p_imag_B, delta_x, delta_y)
    phase_A = ts.get_wave_function_phase(p_real_A, p_imag_A)
    phase_B = ts.get_wave_function_phase(p_real_B, p_imag_B)
    
    return [tot_energy, norm_A, norm_B, phase_A[dim/2, dim/2], phase_B[dim/2, dim/2]]

    
# this is the basic evolution step
def evolution_step():
    ts.evolve_2GPE(p_real_A, p_imag_A, p_real_B, p_imag_B, particle_mass_A, particle_mass_B, external_pot_A, external_pot_B,
                   omega, rot_coord_x, rot_coord_y, delta_x, delta_y, delta_t, iterations,
                   [g_A, g_B, g_AB, Omega_Rabi_real, Omega_Rabi_imag], kernel_type,
                   periods, imag_time)
                   
def calculate_mean_ext_pot():
    sum_norm_A =0.
    sum_norm_B =0.
    sum_pot_A = 0.
    sum_pot_B = 0.
    for yi in range(0, dim):
        for xi in range(0, dim):
            sum_norm_A += p_real_A[xi,yi]*p_real_A[xi,yi]+p_imag_A[xi,yi]*p_imag_A[xi,yi]
            sum_norm_B += p_real_B[xi,yi]*p_real_B[xi,yi]+p_imag_B[xi,yi]*p_imag_B[xi,yi]
            sum_pot_A += (p_real_A[xi,yi]*p_real_A[xi,yi]+p_imag_A[xi,yi]*p_imag_A[xi,yi])*external_pot_A[xi,yi]
            sum_pot_B += (p_real_B[xi,yi]*p_real_B[xi,yi]+p_imag_B[xi,yi]*p_imag_B[xi,yi])*external_pot_B[xi,yi]
        
    return sum_pot_B / sum_norm_B - sum_pot_A / sum_norm_A

Initialize the state

In [35]:
# width of the gaussian envelope
width = 4

def gauss_state(x,y):
    val = np.exp(-(x**2 + y**2)/(2. * width**2))
    return val

p_real_A = np.zeros((dim,dim))
p_imag_A = np.zeros((dim,dim))
p_real_B = np.zeros((dim,dim))
p_imag_B = np.zeros((dim,dim))

for i in range(0, dim):
    y = y_vec[i]
    for j in range(0, dim):
        x = x_vec[j]
        p_real_A[i, j] = np.real(gauss_state(x,y))
        p_imag_A[i, j] = np.imag(gauss_state(x,y))
        p_real_B[i, j] = np.real(gauss_state(x,y))
        p_imag_B[i, j] = np.imag(gauss_state(x,y))

# normalize the initial condition
norm=np.sqrt(delta_x * delta_y * np.sum(np.abs(p_real_A)**2 + np.abs(p_imag_A)**2))
p_real_A=p_real_A/norm
p_imag_A=p_imag_A/norm
norm=np.sqrt(delta_x * delta_y * np.sum(np.abs(p_real_B)**2 + np.abs(p_imag_B)**2))
p_real_B=p_real_B/norm
p_imag_B=p_imag_B/norm

Find the ground state

In [36]:
omega = 0.
imag_time = True    # True: imaginary time evolution; False: real time evolution
delta_t = 0.5e-4     #evolution time of a single iteration step
iterations = 5000   #number of iterations
start_it=1
max_it = 5

g_A =coupling_const
g_B = coupling_const
g_AB = 0.99 * g_A


Omega_Rabi_real = 0
Omega_Rabi_imag = 0


# plot_evolution variables
Time = np.zeros(max_it - start_it + 2)
Tot_energy = np.zeros(max_it - start_it + 2)
Norm_A = np.zeros(max_it - start_it + 2)
Norm_B = np.zeros(max_it - start_it + 2)



time = 0.
[tot_energy, norm_A, norm_B, a, b] = get_info(time)
Time[0] = time
Tot_energy[0] = tot_energy
Norm_A[0] = norm_A
Norm_B[0] = norm_B

for cont in range(start_it-1, max_it):
    evolution_step()
    time=(cont + 1) * iterations * delta_t
    [tot_energy, norm_A, norm_B, a, b] = get_info(time)
    Time[cont+1] = time
    Tot_energy[cont+1] = tot_energy
    Norm_A[cont+1] = norm_A
    Norm_B[cont+1] = norm_B

# calculate overlap between p_A and p_B (only for real p_A and p_B); this is used in the post-processing
k = 0.
for i in range(0, dim):
    for j in range(0, dim):
        k += p_real_A[i,j] * p_real_B[i,j]

k *= delta_x * delta_y
plot_evolution()
/usr/local/lib/python2.7/dist-packages/ipykernel/__main__.py:69: DeprecationWarning:

using a non-integer number instead of an integer will result in an error in the future

Real time evolution

In [37]:
imag_time = False    # True: imaginary time evolution; False: real time evolution
iterations = 200   #number of iterations
start_it=1
max_it = 1600



# set Rabi coupling
Rabi_period = 20. 
Omega = 2. * np.pi / Rabi_period 
delta = np.pi / 3.
Omega_Rabi_imag = 2. * np.pi/Rabi_period * np.sin(delta)
Omega_Rabi_real= 2. * np.pi/Rabi_period * np.cos(delta)

# vectors to be post-processed
Time = np.zeros(max_it - start_it + 2)
Tot_energy = np.zeros(max_it - start_it + 2)
Norm_A = np.zeros(max_it - start_it + 2)
Norm_B = np.zeros(max_it - start_it + 2)
Phase_A = np.zeros(max_it - start_it + 2)
Phase_B = np.zeros(max_it - start_it + 2)
mu_diff_pot = np.zeros(max_it - start_it + 2)

time = 0.
[Tot_energy[0], Norm_A[0], Norm_B[0], Phase_A[0], Phase_B[0]] = get_info(time)
Time[0] = time
mu_diff_pot[0] = calculate_mean_ext_pot()

for cont in range(start_it-1, max_it):
    if cont%100 == 0:
        print cont
    evolution_step()
    time=(cont + 1) * iterations * delta_t
    [Tot_energy[cont+1], Norm_A[cont+1], Norm_B[cont+1], Phase_A[cont+1], Phase_B[cont+1]] = get_info(time)
    Time[cont+1] = time
    mu_diff_pot[cont+1] = calculate_mean_ext_pot()

plot_evolution()
0
100
200
300
400
500
600
700
800
900
1000
1100
1200
1300
1400
1500
/usr/local/lib/python2.7/dist-packages/ipykernel/__main__.py:69: DeprecationWarning:

using a non-integer number instead of an integer will result in an error in the future

Post processing

In [38]:
print mu_diff_pot
[ -1.77635684e-15   5.29765032e-09   4.23889635e-08 ...,   6.89187752e-01
   6.91339276e-01   6.93516060e-01]

Calculate $\eta = \frac{N_2 - N_1}{N}$, $\phi_{21} = \phi_2 - \phi_1$, $\frac{d}{dt} \eta$ and $\frac{d}{dt} \phi_{21}$

In [49]:
eta = np.zeros(max_it - start_it + 2) 
dt_eta = np.zeros(max_it - start_it + 2) 
phase_21 = np.zeros(max_it - start_it + 2) 
dt_phase_21 = np.zeros(max_it - start_it + 2) 

eta  = 0.5 * (Norm_B - Norm_A)
phase_21 = Phase_B - Phase_A

# calculate time derivative of eta and phi_21
for i in range(1,max_it - start_it + 1):
    dt_eta[i] = (eta[i+1] - eta[i-1]) / (2 * delta_t * iterations)
    var = phase_21[i+1] - phase_21[i-1]
    if np.sign(var) * (var) > np.pi:
        var = var - np.sign(var)*2*np.pi
    dt_phase_21[i] = (var) / (2 * delta_t * iterations) 

dt_eta[0]=dt_eta[1]
mu_1 = g_A / length**2 + g_AB / length**2
mu_2 = g_B / length**2 + g_AB / length**2

_dt_eta = -2. * Omega * k * np.sqrt(1. - eta**2) * np.sin(delta + phase_21)   # derivative of eta deduced from eta and phi_21
_dt_phase_21 = - (g_A-g_AB) / length**2 * eta + 2 * Omega * k * eta / np.sqrt(1. - eta**2) * np.cos(delta + phase_21)

Plot $\frac{d}{dt}\eta$ (with plotly)

In [50]:
x = np.arange(0, 1001 * delta_t * iterations, delta_t * iterations)
trace1 = go.Scatter(
    x = x,
    y = dt_eta,
    mode = 'lines',
    name = 'eta from simulation'
)

trace0 = go.Scatter(
    x = x,
    y = _dt_eta,
    mode = 'lines',
    name = 'eta from equation'
)
data = [trace0, trace1]

layout = dict(
    font = go.Font(size=14),
    xaxis=dict(
        title='Time',
        showticklabels=True,
        autotick=False,
        ticks='outside',
        tick0=0,
        dtick=0.5,
        ticklen=10,
        #tickwidth=4,
    ),
    yaxis = dict(
        title = '',
        autotick=False,
        ticks='outside',
        tick0=0,
        dtick=0.25,
        ticklen=10,
        #tickwidth=4,
    ),
)

fig = dict(data=data, layout=layout)

iplot(fig)
#py.image.save_as(fig, 'eta-time-derivative-null-deltapi3.pdf') 
Drawing...

Plot $\frac{d}{dt}\phi_{21}$ (with plotly)

In [51]:
x = np.arange(delta_t * iterations, max_it * delta_t * iterations, delta_t * iterations)
trace1 = go.Scatter(
    x = x,
    y = dt_phase_21,
    mode = 'lines',
    name = 'phi_21 from simulation'
)

trace0 = go.Scatter(
    x = x,
    y = _dt_phase_21,
    mode = 'lines',
    name = 'phi_21 form equation'
)

data = [trace0, trace1]

layout = dict(
    font = go.Font(size=14),
    xaxis=dict(
        title='Time',
        showticklabels=True,
        autotick=False,
        ticks='outside',
        tick0=0,
        dtick=0.5,
        ticklen=10,
        #tickwidth=4,
    ),
    yaxis = dict(
        title = '',
        autotick=False,
        ticks='outside',
        tick0=0,
        dtick=0.25,
        ticklen=10,
        #tickwidth=4,
    ),
)

fig = dict(data=data, layout=layout)


iplot(fig)
#py.image.save_as(fig, 'phi_21-time-derivative-null-deltapi3.pdf')
Drawing...
In [52]:
print - (g_A-g_AB)/(length*length)*eta
[  9.53989416e-18   8.50116914e-05   1.70020470e-04 ...,  -5.09517865e-04
  -5.92551080e-04  -6.75561184e-04]
In [78]:
print Norm_A[250]
1.87662399815
In [88]:
x = np.arange(0, 1001 * delta_t * iterations, delta_t * iterations) * 0.1
trace0 = go.Scatter(
    x = x,
    y = dt_eta,
    mode = 'lines+markers',
    marker = dict(maxdisplayed = 62),
    name = '$  \eta $'
)

trace1 = go.Scatter(
    x = x,
    y = _dt_eta,
    mode = 'lines',
    name = '$ f_{\dot{\eta}}(\eta, \phi) $'
)

trace2 = go.Scatter(
    x = x,
    y = dt_phase_21,
    mode = 'lines+markers',
    marker = dict(maxdisplayed = 62),
    name = '$  \phi $',
    yaxis='y2'
)

trace3 = go.Scatter(
    x = x,
    y = _dt_phase_21,
    mode = 'lines',
    name = '$ f_{\dot{\phi}}(\eta, \phi) $',
    yaxis='y2'
)

data = [trace0, trace1, trace2, trace3]

layout = dict(
    legend=dict(
        x=0.01,
        y=1
    ),
    width = 900,
    height = 500,
    margin = dict(t=0),
    font = go.Font(size=15),
    xaxis=dict(
        title='$t/T$',
        showticklabels=True,
        autotick=False,
        ticks='outside',
        tick0=0,
        dtick=0.25,
        ticklen=10,
        #tickwidth=4,
    ),
    yaxis = dict(
        title = '$ \eta $',
        autotick=False,
        ticks='outside',
        tick0=0,
        dtick=0.25,
        ticklen=10,
        #tickwidth=4,
    ),
    yaxis2 = dict(
        title = '$  \phi $',
        autotick=False,
        ticks='outside',
        tick0=0,
        dtick=0.25,
        ticklen=10,
        overlaying='y',
        side='right',
    ),
)

fig = dict(data=data, layout=layout)

iplot(fig)
#py.image.save_as(fig, 'phi.pdf')
Drawing...
In [ ]: